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Abstract 

We describe in detail two numerical simulation methods valid to study systems 
whose thermostatistics is described by generalized entropies, such as Tsallis. The 
methods are useful for applications to non-trivial interacting systems with a large 
number of degrees of freedom, and both short-range and long-range interactions. 
The first method is quite general and it is based on the numerical evaluation of the 
density of states with a given energy. The second method is more specific for Tsallis 
thermostatistics and it is based on a standard Monte Carlo Metropolis algorithm 
along with a numerical integration procedure. We show here that both methods 
are robust and efficient. We present results of the application of the methods to 
the one-dimensional Ising model both in a short-range case and in a long-range 
(non-extensive) case. We show that the thermodynamic potentials for different val- 
ues of the system size N and different values of the non-extensivity parameter q 
can be described by scaling relations which are an extension of the ones holding 
for the Boltzmann-Gibbs statistics (q = 1). Finally, we discuss the differences in 
using standard or non-standard mean value definitions in the Tsallis thermostatis- 
tics formalism and present a micro canonical ensemble calculation approach of the 
averages. 
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1 Introduction 



Non-extensive systems are those for which the thermodynamic potentials do 
not scale linearly with the system size. As a way of example, in some electric or 
magnetic systems with very long-range interactions the ground state energy 
per particle increases with the number of particles. In the absence of other 
effects, such as screening, these system are "genuinely" non-extensive. If we 
apply to them the standard Boltzmann-Gibbs formalism of the Statistical 
Mechanics, we find that the internal energy, Helmholtz free energy and other 
thermodynamic potentials are non-extensive as well. This standard formalism 
can be implemented by using the definition of the entropy S in terms of the 
probabilities pi of the i = 1, . . . ,W possible microscopic configurations^: 

S = -^Pilnpi (1) 



The actual calculation of the entropy assumes a set of probabilities {pi}- These 
are computed by finding the maximum of the above expression when some 
extra conditions defining an ensemble (fixed number of particles and mean 
energy, for example) are imposed. 

Even for systems in which the energy levels do scale with the system size, 
it is possible, by using generalized definitions of the entropy, to obtain non- 
extensive thermodynamic potentials. One of the most successful generaliza- 
tions is that of Tsallis which in 1988 proposed [1] the following alternative 
expression for the entropy: 

S q = ^ ( 2 ) 



where q is an entropic index that characterizes the degree of non-extensivity. 
It is possible to show that the entropy of the composed system A+B satisfies 
the relation: 

S q (A + B) = S q (A) + S q (B) + (1 - q)S q (A)S q (B) (3) 



when A and B are independent systems in the sense that Pij(A + B) = 
Pi(A)pj(B). We see that for q ^ 1 there is no additivity in the entropy, which 
also implies non-extensivity. The Boltzmann-Gibbs entropy, Eq. (1), and ex- 
tensivity are recovered in the limit q — > 1. Since the probabilities {pi} satisfy 
p\ > pi for q < 1 and pf < pi for q > 1, the superextensive, q < 1, and 

4 We use throughout the paper dimensionless units where the Boltzmann constant, 
ks is equal to 1 
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the subextensive, q > 1, regimes will privilege the rare and frequent events 
respectively. 

In the last years there have been many studies in which Tsallis non-extensive 
statistics has been applied to different situations (See [2,3] for a review). 
In some cases, the systems considered are genuinely non-extensive (in the 
sense defined above) while in others the non-extensivity arises as a result of 
the application of the new statistics. In fact, and due to the intrinsic non- 
extensivity of the Tsallis statistics, it has been argued that its natural range 
of applicability should include systems with long-range interactions or long- 
range microscopic memory processes, as well as dynamical systems in which 
the space-time geometry has a mult ifr act al-like structure, because those sys- 
tems are in general genuinely non-extensive. Although most of the literature 
(including this paper) basically derives equilibrium properties starting from 
the generalized definition of the entropy, it has been conjectured recently [3], 
however, that the Tsallis entropy could be relevant instead in the study of 
non-equilibrium processes. 

Due to the difficulty of deriving exact results, it is natural to use numerical 
methods to obtain the properties of a system with many degrees of freedom 
when studied under the rules of the new statistics. This is necessary in order 
to extract results that could be checked against experiments. However, these 
studies have been hampered by the failure of the typical Monte Carlo methods 
to adequately generate representative equilibrium configurations distributed 
according to Tsallis statistics. It is the purpose of this paper to explain in detail 
new methods that can be used to study the equilibrium properties of a many- 
particle system when it is considered under generalized statistics. Although 
our methods are quite general, we will illustrate their use by considering a pro- 
totypical genuinely non-extensive system: the Ising ferromagnet model with 
long-range interactions. We will also consider the short-range Ising model in 
order to test the simulation methods and to compare the results obtained from 
the use of Tsallis statistics in extensive and non-extensive systems. 

In the remaining of the section, we will outline briefly which are the basic dif- 
ficulties one encounters when trying to generalize the standard Monte Carlo 
methods (such as the Metropolis algorithm) to the study of Tsallis statistics. 
The main problem is that the probabilities {p^ can not be given an explicit 
expression, as we will see in the following discussion. Let us consider the canon- 
ical ensemble. The probabilities {p^ in this ensemble are found by solving the 
maximization problem of the entropy S q as given by Eq.(2) subject to the 
constrains of (i) positivity: pi > 0, (ii) normalization: J2iPi — 1 an d (iii) a 
fixed mean value for the internal energy: (Ti) = U, where 7i is the Hamilto- 
nian of the system and the mean value of any function O of the microscopic 
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configurations is computed according to the general rule: 



(0> = SOi«(p i ), (4) 



Oi is the value of O at the configuration whose probability is pi and we have 
introduced a function u(pi) that allows the definition of generalized mean val- 
ues. The standard mean values are recovered by taking u(pi) = pi. Although, 
initially the choices u{j>i) = pi {first option), and u(pi) = pf (second option) 
were considered, later, it was shown that a better choice, in the sense that it 
preserves the Legendre structure of the resulting thermodynamics formalism, 
is to consider u(pi) = Pi/Y^jPj (third option) [4]. We will use the following 
notation of the averages in this third option: 

(0) ? = E^; p i = ^~ q - (5) 



{Pi} are known as the "escorts" probabilities [5]. It is possible to recover the 
configuration probabilities pi from the escorts probabilities using: 

f. = ^rpr.- («) 

^3 3 



The entropy, in terms of the {-Pjj's is given by: 

1 _ (y.p 1/q )~i 

s q = ^ ; ' (7) 



Concerning the different definitions for the averages, it should be said that it 
has been shown recently [6-8] that the standard mean values of the first op- 
tion, u(pi) = Pi, can be also made compatible with the Legendre structure of 
the thermodynamics and the resulting formalism also represents a thermody- 
namically stable description. In this paper, we follow mainly the formulation 
in terms of the mean values defined by Eq.(5), although in a later section we 
will show that the results obtained using the standard mean values can be 
mapped onto the ones obtained using Eq. (5). 

The maximization problem for the unknown escort probabilities Pi in the 
canonical ensemble with a given internal energy U q is: 



L-[S q -f3j2e t P l -aJ2P l ]=0 (8) 
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Pi>0 

i 



(9) 
(10) 

(11) 



where a, f3 are Lagrange multipliers. {ei\ are the energy levels of the system 
under consideration whose ground state energy will be denoted by Eq. Solving 
the problem Eq.(8-ll) one obtains the probabilities for the canonical ensemble 
as [4]: 



P = < 



0, 



1 



[i-a-9M^-u q )/(E j pj /9 ) q ] TJ - 



(l-q)f3(Ei-U q ) 

(E 3 ^ /9 ) 9 



< 



otherwise 



(12) 



The probabilities defined in this way are real and non-negative. The condition 
giving the possible values of j3 and E{ for which Pi ^ in Eq. (12) is called the 
cut-off condition. One can show that the probabilities Eq. (12) are invariant 
under a change in the energy levels £j — > + Ae (and the same change in 
the internal energy U q — ^ U q + As) for arbitrary Ae. By introducing the 
temperature T — it is possible to show also the validity of the relation 
[7,4] 

!/T = dS q /dU q (13) 



which reflects the Legendre structure of the thermodynamics obtained. 

Notice that Eq.(12) does not give yet the actual values of the probabilities 
since the {-Pjj-'s appear in a non-trivial way in both sides of the equation 
(either explicitly or in the cut-off condition). This is different from the solution 
obtained in the usual Boltzmann-Gibbs canonical ensemble (recovered in the 
limit q — > 1) in which the solution adopts the explicit form: 

P = ^— w ( 14 ) 

(although, of course, it is very difficult to compute the denominator of this 
expression, the partition function, for an interacting system). An iterative 
method to solve Eq.(12) has been used in [4]. In this method, an initial guess 
for the probabilities is fed in the right hand side of (12) and this equation is 
used recursively until convergence is achieved. We will see, however, that for 
many-particle systems, it might very difficult to achieve convergence in some 
cases. 
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A convenient way of writing Eq. (12) is by using an auxiliary parameter (5' 
defined as [4]: 



(l-q)PX j e j P j + G: j P} /q )-* 



P' = ~ — g ~ - ( 15 ) 



Defining T' = and using Eqs. (7) and (11) one can rewrite this equation 



as 



T = r-v-<)v, (16) 

1 + (1-9)S, 1 ' 



In terms of (5' the solution adopts a form similar to that of the standard 
canonical ensemble: 



P = { 



0, 1 - (1 - q)P'e i < 



H-d-^e^ otherwise 



One can then adopt the following practical procedure [4]: choose a value for 
the parameter T' and compute the probabilities Pi as a function of T' us- 
ing Eq.(17). Compute the internal energy U q , the entropy S q and the tem- 
perature T, always as a function of T', using Eqs. (11), (7) and (16), respec- 
tively. Finally, vary T' in order to make the parametric plots U q (T) and S q {T). 
Other thermodynamic potentials follow the usual definition. For instance, the 
Helmholtz free energy is F q = U q — TS q . 

It is important to realize that although the probabilities Pj, when considered 
as a function of T, do not depend on an arbitrary shift Ae of the energy levels 
or, in other words, do not depend on the zero of energy, Eq, they do depend 
on Eq when considered as a function of T'. This means that the averages as a 
function of T' can not be physically relevant because they depend on the zero 
of energy. This is why T' has to be interpreted only as an auxiliary parameter, 
not as an actual temperature. Of course, the relation V — > T depends also 
on the zero of energy in such a way that both dependences cancel and the 
averages as a function of T are independent of a shift in the energy levels. An 
interesting question is to determine the range of values for the parameter T' 
that should be used in order to obtain the usual range < T < oo for the 
temperature T. Using that, according to the definition (2), it is 1+(1— q)S q > 
we obtain from Eq.(16) that V should vary in the range [(1 — q)E , oo), where 
we have used that the energy at zero temperature is the ground state energy 
U q (T = 0) = Eq. Therefore, it is important to use the right range of values 
for T' in order to reach all the possible values for T. In particular, T' might 
need to take negative values either for q < 1 or for q > 1 unless one adopts 
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Eq = as we will throughout this paper. To our understanding, it is not clear 
in the literature the fact that the averages as a function of the parameter 
T' depend on the zero of energy and that it might be necessary to consider 
negative values for T' in order to span the whole range of values for T. 

As stated before, the main problem to perform Tsallis thermostatistics simu- 
lations at a given temperature T is that there is not an explicit expression for 
the probabilities Pj, c.f. see Eq. (12). The practical procedure outlined above 
(compute Pi as a function of T' and then compute U q , S q and T as a function 
of T' in order to make parametric plots by varying T') is not straightforward 
to implement numerically since it is very difficult to use Eq.(7) to compute the 
entropy. This is because the usual Monte Carlo methods, i.e. the Metropolis 
algorithm, require only the probabilities Pi up to a normalization factor where 
Eq.(7) requires the absolute, normalized values of Pj. It is the object of this 
paper to explain in detail some numerical methods of the Monte Carlo type 
that can be used to perform the necessary averages for generalized statistics, 
including Tsallis. 

There have been previous attempts to perform numerical simulations of Tsallis 
statistics using Monte Carlo methods. An earlier work in this direction is that 
of T. Penna et. al. [9] who extended the Metropolis acceptance procedure 
to include a dependence in the q parameter. However, this method does not 
satisfy the detailed balance condition which is a key ingredient of Monte Carlo 
methods. Another interesting approach is that of I. Andricioaei et. al. [10], who 
performed a Metropolis Monte Carlo algorithm which does satisfy the detailed 
balance condition for the probability P, as a function of the parameter T' but, 
since they do not make the temperature transformation T' — > T, they are 
unable to determine the actual temperature T of the simulation. All these 
works considered the second version, u(j>i) = pj, for the definition of averages, 
which, as discussed before, has proven afterwards not to be the optimal election 
[4]. We have used also a similar sampling in the context of Simulated Annealing 
[11,12]. A recent approach proposed by A. R. Lima et. al. [13] uses the broad 
histogram Monte Carlo method, which determines the number of microstates 
using a balance equation between near neighbor energy levels. They are able 
to apply this method to the Ising Model with short-range interactions. This 
is a valid Monte Carlo simulation with full control of the temperature T but 
its applicability is somewhat restricted. As we will show, the Ising model with 
long-range interactions can not be treated straightforwardly with this method 
because the spin flip dynamics does not produce transitions between near 
energy levels and, consequently, the broad histogram Monte Carlo method, in 
its present form, can not be used to study long-range interacting systems. We 
are able to overcome these difficulties by extending a method which had been 
developed some time ago [14-17] to compute directly the number of states with 
a given energy and which does not depend on the definition of the entropy. In 
fact, our method can be used to study any statistics and applications will be 
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shown both for the Boltzmann-Gibbs and the Tsallis statistics. We develop 
yet a second method which is devised specifically for the Tsallis statistics and 
that has the advantage that it uses the familiar Metropolis algorithm plus a 
numerical integration. 

This paper is organized in the following form: in section 2 we use a simple and 
limited enumeration procedure valid only for small system sizes. However, this 
method is exact and can be used to check against some of the approximated 
methods we will introduce later. In this section we also introduce the short- 
range Ising model (SRIM) and the long-range Ising Model (LRIM). These two 
models will be employed in this paper in order to test the numerical meth- 
ods described here and to compare the behavior of the non-extensive Tsallis 
thermostatistics in genuinely extensive and non-extensive systems. In section 
3 we explain in some detail the Histogram by Overlapping Windows (HOW) 
method. We show that this method has a wide range of applicability since it 
can be used both for short-range and long-range systems as well as for any 
kind of statistics. Section 4 presents a Metropolis Monte Carlo type method 
specially devised to study systems in the Tsallis statistics. In section 5 we 
present some results concerning the validity of some scaling relations for the 
SRIM and LRIM in the Tsallis thermostatistics context. These scaling rela- 
tions are an extension of the ones holding for the Boltzmann-Gibbs statistics 
. In section 6 we present results for these models using standard mean values 
instead of those defined by Eq. (5). In section 7 we discuss the results of us- 
ing Tsallis statistics in the microcanonical ensemble. Finally, in section 8 we 
summarize the main conclusions of this work. 



2 Exact enumeration 



The problem to compute the probabilities {Rj} using Eq.(17) is that the num- 
ber of terms in the sum of the denominator of this equation, the number of 
microstate configurations W, is extremely large (typically scales exponentially 
with the system size). However, for small systems, it might be possible to enu- 
merate completely the microstates and, therefore, to compute magnitudes of 
interest such as the internal energy U q (T), the entropy S q (T), the free energy 
F q (T), etc. We follow this approach in this section. Although the system sizes 
one can usually study with this method are very far from reaching a situation 
in which scaling laws with system size apply, we use the results as a bench-test 
in order to compare with other approximate methods that will be introduced 
in the following sections. 
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We will consider Ising type models with Hamiltonian: 

n = Y,Mi-*i*i) ( 18 ) 

(i,3) 

where each of the A" spin variables Sj, % — 1, . . . , TV can take the values ±1. The 
sum J2(i,j) runs over all distinct pairs of sites on a (/-dimensional regular lattice 
of linear size L = N l l d with periodic boundary conditions and lattice constant 
equal to 1. is the coupling parameter between spins % and j. Note that for 
ferromagnetic couplings, > 0, the ground state is double degenerate and 
its energy is E = 0. The usual, nearest-neighbors, or short-range Ising model 
(SRIM), is obtained taking = 1, if = 1 and = 0, if > 1. The 
long-range Ising model (LRIM) is defined by using 

Jy = 1/r?-, (19) 



where is the distance between the spins i and j, and the parameter a sets 
the interaction range. The SRIM is formally recovered by taking the limit 
a — > oo. Depending on the value of a and the space dimension d, the LRIM 
has two regimes: the extensive regime, a > d, and the non-extensive regime, 
a < d. This can be seen by roughly estimating the mean energy per spin in 
an infinite system as f£° drr d ~ 1 r~ a . We obtain a convergent integral if a > d 
(extensive behavior), and for a < d the integral diverges (non-extensive). 
More precisely, a convenient scale for the mean energy per spin in a finite 
system of size N is given by [18]: 

N = 1 + d [ drr^r-" = ^H^lA (20 ) 
J 1 — a/d 



The definition of N is such that the limit a — > d is well defined. Again, we see 
that for a > d the internal energy per spin scales as a constant in the limit 
of large N, but for a < ad, it grows with the system size. The system is, in 
this latter case, genuinely non-extensive. The SRIM limit, a — > oo, gives the 
expected result N — 1. 

The number of configurations in the Ising model is W = 2 N . We have made 
a complete enumeration of the % — 1, . . . , W configurations and their energies 
£j for a linear chain, d — 1, of sizes up to = 34. We have used these results 
to compute the probabilities P« and then the thermodynamic magnitudes of 
interest using Tsallis statistics with the third option for the averages. In Fig. 
1, we plot the exact internal energy U q as a function of the temperature T, for 
several values of the parameter q for the LRIM in a genuinely non-extensive 
situation, a = 0.8, Fig. l.b, and for the SRIM, Fig. l.a. In the q < 1 case we 
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observe that there is a range of temperatures for which the internal energy is 
not a single-valued function of the temperature and, for a given value of T, 
there are several possible values for U q . This ambiguity is resolved by using 
a Maxwell-like construction [19,4] that replaces the loop in the energy curves 
by a vertical straight line connecting the two points with the same free energy 
F q (T). 




q=0.6 
q=1.0 
q=1.4 



q=0.6 
q=1.0 
q=1.4 



-20 -10 10 20 

logT 

Fig. 1. Internal energy U q as a function of the temperature T, from the exact 
evaluation of Eqs. (11,12) for one-dimensional Ising models with N = 34 spins. 
Plot (a) is for the SRIM and plot (b) for the LRIM in a genuinely long-range case, 
a = 0.8. The insert shows the free energy F q for q = 0.6. In this case, a part of the 
energy curve is replaced by a vertical straight line (shown by dots in the plot) such 
that the Maxwell criterion of equal free energies is satisfied. 

We stress that the loops in the energy curves appear as a result of the T' — > T 
transformation and, therefore, will not be observed when plotting the energy 
as a function of T' . Typical T' — > T transformations are shown in Fig. 2. We 
observe in this figure that for q < la same value of T' can correspond, in some 
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-10 10 20 

log T' 




-10 10 20 

log r 



Fig. 2. The T' — > T transformation for the one-dimensional Ising models as obtained 
from the exact evaluation of Eq. (16) for a system size N = 34 and different values of 
q. Plot (a) is for the SRIM, while plot (b) is for the LRIM for a = 0.8. Notice that, in 
both cases, there is a temperature range for which the curves are almost horizontal. 
This makes it very difficult to use iterative methods for the determination of the 
probabilities using directly Eq. (12). 

cases, to three or more values of T giving rise to the observed multivalued 
behavior in the energy. Figure 2 helps us to understand the failure of the 
iterative method that has been proposed [4] to solve the set of equations (12). 
Although each value of V defines a unique set {Pi}. We see Fig. 2 that for 
q > 1 there are some intervals of T' where the transformation T' — > T is 
almost horizontal. Therefore, one value of T corresponds nearly to a complete 
interval for T 1 and hence, there are many possible solutions for {Pi} very close 
to the real one. This is the main reason for the failure of iterative methods for 
q > 1. The situation worsens for increasing system size N. 
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Whatever illuminating the method of exact enumeration is, its validity is lim- 
ited to very small values for N. To our knowledge, the largest value ever 
considered in an exact enumeration scheme for a short-range Ising model is 
N = 4 3 = 64 [20]. Simulations at larger N sizes require other methods, as the 
ones presented in the next two sections. 



3 The Energy Histogram Method using Overlapping Windows 



Although the number of possible microscopic configurations is in general very 
large, the range of possible energy values usually takes a much smaller value. 
For instance, for the one-dimensional SRIM introduced in the previous section, 
with N spins we have W — 2 N , but the number of possible energy values is 
N/2 + 1. Let M, in general, be the number of possible energy levels. We 
denote by fl(E k ) the number of microscopic states sharing the same energy 
Ek, k = 0, . . . , M — 1. Obviously J2k &(Ek) = W. We rewrite all the sums in 
Eqs. (17,7,5) as: 



0, 1 - (1 - q)!3'E k < 
-^fcl 1 ^ g _ ? ot herwise 

-£ n n(E n )[l-(l-q)i3'E n ]T^ 

s = l-EtfllW) 1 ^)"' (22) 

(0) q = J2tt(E k )0(E k )P(E k ) (23) 

k 



where the sums run over the M energy levels. 

Notice that, once the f2(£^)'s have been computed, any statistics can be per- 
formed upon the system. Whether we use Tsallis, Boltzmann-Gibbs or any 
other generalized statistics is simply a trivial change in the computational 
scheme. Moreover, it is also trivial to compute the averages for any value of 
the parameters, say T or q. Therefore, although the calculation of the fl k s 
might be time consuming, the pay-off is tremendous^. 

In general, the Q(E k ) are very difficult to obtain exactly. An important excep- 
tion that will be used throughout this paper is the SRIM in 1-d, for which the 
energy levels are given by E k = Ak for k — 0, . . . , N/2 and whose degeneracy 



5 All the simulations in this paper have been performed using a Pentium-Ill pro- 
cessor at 550MHz 
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is: 



(24) 



In the cases in which Q(E k ) is not known we need approximate numerical 
methods. The most naive way to find the f^-E^'s is to generate different 
system configurations randomly and count how many times a configuration 
with energy E k appears. However, this approach fails because the complete 
set of fl(Ek) values span too many orders of magnitude. In general, two energy 
levels E k and E n could differ as much as fl k /fl n ~ expiV. This means that 
the range of variation of Vt(E k ) over the M different energy levels is very large 
and it is not possible to generate in a single run a histogram that covers all 
the energy levels, unless one generates a number of configurations of the order 
of the total number available to the system, W . 

The Histogram by Overlapping Windows method (HOW) used here [14] avoids 
this problem by generating system configurations within a restricted energy 
interval and estimating the relative weights of these energy levels from the 
number of times they appear in each interval. By generating enough intervals 
spanning the whole energy range, one is able to obtain good quality estimators 
of the numbers Q(E k ). An earlier account of the method has been given in 
[21] and we explain now in some detail how the method works. 

Let us consider first the SRIM in arbitrary dimension. In this case, the possible 
energy values are E k = 4k for k — 0, . . . , dN/2. Following the original work 
[14], we consider the intervals (windows) {Eq, E\, E2, E%}, {£3, E4, E$, Eq}, 
{Eq, E7, Eg, Eg}, etc. Each window consists of 4 consecutive energy levels and 
the last energy value of one window is the first of the next one. The next step 
is to take one of the intervals and to generate configurations whose energy 
belongs to it. This is achieved, after preparing the system initially with one of 
the energies of the interval, by flipping spins chosen at random. A spin flip is 
accepted only if it leaves the system in one of the energy levels of the interval 
and it is rejected otherwise. The ratio of the number of generated states with 
energy E k to the number of generated states with energy E n is an unbiased 
estimator of Vt(E k ) /Q(E n ), for E k and E n within the energy window. The 
quality of the estimator increases with the number of generated configurations. 
From the overlap between windows one can compute fl(E k ) for the whole range 
of energies. The number of energy values in each window (4 in the previous 
example) is not important as far as it is not too large (such that the ratios 
Q(E k )/Q(E n ) are not extremely small or large) and it is not too small either. If 
the window is very small, most spin flips will yield an energy outside the range 
of allowed values and the number of accepted, i.e. independent, configurations 
will be very small. Moreover, the final algorithm must be ergodic: any energy 
value in a window should be obtained from any other value in the same window 
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after a sufficient number of spin flips. 
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Fig. 3. The number of states £l(Ek) for N = 34 calculated exactly (lines), and using 
the HOW method (symbols), (a) SRIM, (b) LRIM. 

The same basic idea has been used in other short-range Hamiltonians [14- 
17,22]. We are concerned now with the extension of this method to consider 
long-range interactions, in particular the LRIM introduced before. A modi- 
fication needed is that the energy values E k will represent now a continuum 
set of energies with a bin size SE. The energy levels are then Ek = kSE 
k — 1, . . . , M and Q(Ek) counts all the configurations % whose energy sat- 
isfies Ek < Si < Ek + SE. In other words, one makes the approximation of 
considering that all the energies lying in one bin count as one single level. This 
turns out not to be a critical point, although one has to check that the results 
are independent, within the simulation errors, of the magnitude of SE. 

A more important point concerns the optimal size / of the energy window 
{Ek, Ek+i, . . . , Ek+i}. Since, for a long-range interacting system, a single spin 
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Fig. 4. All the possible jumps from the energy level mark in black generated by a 
single spin flip, (a) SRIM, (b) LRIM. Note in (b) the small modulation in 0,(E^) of 
the order of l E = 2J2iLi r~ a (See the text). 



flip can produce a very large change in the energy, it is important not to choose 
/ too small. To make this point clear, we have calculated exactly the number of 
states Q(Ek) for a system with N = 34 by using the complete enumeration of 
the W = 2 34 possibles configurations, see Fig. 3. Using these exact results we 
study the energy changes that a single spin flip makes both in the SRIM and 
the LRIM cases. A typical situation is shown in Fig.4. In this figure we plot 
the histogram of the (exact) number of configurations using a value for the 
bin SE = 4 for the SRIM, and SE = 1 for the LRIM. We select several config- 
urations belonging to one of the energy bins (marked black in the figure) and 
we dash all the levels that are obtained from these configurations by nipping 
one spin. We see that, as expected, the change in energy for the SRIM brings 
the system to one of its neighboring energy levels. However, for the LRIM 
the energy changes are very large and, in fact, the nearest-neighbor energy 
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Fig. 5. The number of states £l{E k ) calculated using the HOW method for the 1-d 
Ising models for several values of N. (a) SRIM, (b) LRIM. In (b) we have shifted 
vertically the curves to avoid overlapping. The inserts show the collapse of all the 
curves using Eq. (25). 

levels can not be reached by using the spin flip dynamics [23]. A measure of 
the typical change in energy obtained when flipping one spin is estimated by 
considering the ferromagnetic ground state configuration with all the spins 
pointing in the same direction. One spin flip in this configuration produces 
a change in energy AE = 2J2iLi r ij a - The equivalent number of energy bins 
is £ = AE/SE. We finally take the size of the energy windows / = 3£. In 
order to make sure that ergodicity is satisfied, we adopt a large overlap of 
size 2£ between the windows. This means that a window goes from E k to 
but the next window goes from E k+ ^ to E k+4 ^ and so on. For the win- 
dow {E k , E k+ i, . . . , E k+3 ^} only those values in the interval {E k+ ^, . . . , E k+2 ^} 
are considered for the evaluation of the ratios Q(E k )/Q(E n ). To summarize, 
for the LRIM, it is necessary that both the window size and the overlap be- 
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tween windows has the correct size, depending on a. For a = 0.8, used in our 
simulations, we take £ = 4 independently of the system size and adjust SE 
accordingly. We have checked that £ = 10 gives the same results within the 
numerical errors. The number of configurations necessary increases with the 
required accuracy. We have adopted in our simulations the criterion that the 
minimum number of counts for any energy bin within a window is 100. The 
knowledge of the exact degeneration of the ground state Q(Eo) = 2 allows 
finally the calculation of all the Q(Ek). 
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Fig. 6. Internal energy U q (T) for 1-dimensional Ising models with N = 34 and 
q = 0.8, 1.0, 1.2. Symbols are obtained using the HOW method and lines show the 
exact results, (a) SRIM, (b) LRIM. 

In Fig. 3, we plot the number of states computed either exactly or by using the 
HOW method, for iV = 34 both for the SRIM and the LRIM. At the resolution 
of the figure, the exact results and the approximate ones are indistinguishable. 
This serves as a test that the HOW method, as implemented here, is capable 
or reproducing accurately the number of states in a known case. In Figs. 5, we 
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Fig. 7. Internal energy U q as a function of temperature T for one-dimensional Ising 
models and different values of N and q, as coming from the application of the HOW 
method in the (a) SRIM and (b) LRIM. The values of q are q = 0.8, 0.9, 1.0, 1.1, 1.2 
(curves from left to right). 

show the number of states for sizes N = 34, 100,200,400, 1000 as computed 
using the HOW method. The inserts of these figures show that this function 
scales as: 



n(E k ) = e N ^/ N ^ 



(25) 



This is valid both for the SRIM and for the LRIM, if we recall that N — 1 for 
the SRIM. The scaling function <j>(x) is different for the SRIM and the LRIM. 
For the SRIM, the result (24) leads to: 

<t>SRiM{x) = X - In (J - l) - In (l - |) (26) 
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No analytical expression is available for the LRIM. Fig. 6 compares the result 
for the internal energy U q obtained using the HOW method and the exact 
results obtained from the exact enumeration algorithm in the case N = 34. 
Again, we can see that differences are too small to show up in this plot. Finally, 
in Fig. 7 we make a similar plot for larger values of N obtained in this case 
by application of the HOW method. 



4 The Monte Carlo method 



We have mentioned already that the usual Monte Carlo algorithms of the 
Metropolis type can not be applied to Tsallis statistics, because the probabil- 
ities {Pi} are not known as a function of the temperature T. However, it is 
possible to use them to compute averages as a function of the parameter T' be- 
cause the probabilities are known as a function of T' except for a normalizing 
factor, which is irrelevant in the Monte Carlo methods. Those averages can 
be performed by using the Metropolis algorithm to generate configurations 
distributed according to the probabilities {Pi} as follows: consider the config- 
uration % with energy flip a spin chosen at random to produce configuration 
j with energy £,-, accept this change with a probability min(l, Pj/ Pi). For the 
Boltzmann-Gibbs statistics the acceptance probability is the celebrated factor 
min[l, exp(— P(ej — £,))]. For Tsallis statistics, Eq. (17), it is instead: 



P(i -> j) 



mm 



, i-(i- g )/y ej - 



9 

1-9 



1 - (1 - q)(3'e j < 
otherwise 



(27) 



This acceptance probability was used for the first time by I. Andricioaei et. 
al. [10]. After generating M configurations, the averages (0) q at fixed T' are 
obtained as the sum: (0) q = A/'- 1 E^iO(s), where O(s) is the value of the 
observable O at the configuration s in the sequence of configurations generated 
by the Monte Carlo method. 

Still, we need to perform the T' — * T transformation, in order to plot averages 
with respect to T, using Eq. (16). In this equation, we can use the Monte Carlo 
averages for the internal energy U q = (Ti) , but the entropy is yet unknown. In 
order to compute the entropy, we combine Eq. (16) and Eq. (13): 

ds, = i + r28 , 

dU q T' — (1 — q)U q 1 ) 
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Fig. 8. The integration procedure implicit in the Eq. (30) to perform the T' — > T 
transformation in the Monte Carlo method. We use here N = 34, q = 0.6 and 
T A = 0. (a) SRIM, (b) LRIM. 

which can be integrated between arbitrary points A and B: 



1 ln[l + (l-q)S l "' 



1-q 



U q (B) 



(29) 



to obtain: 



Uq(B) 



S q (B) = 



[l + (l-q)S q (A)]e 



Uq(A) 



T'-(l-q)Uq 



-1 



1-q 



(30) 



20 




-12 -8 -4 



logT 



160 i 




-12 -8 -4 



logT 

Fig. 9. Internal energy U q as a function of the temperature T, for N = 34 and 
q = 0.6. By symbols we show the results of the Monte Carlo method and by lines 
the exact results. In the inserts we show the free energy criterion explained in the 
main text, (a) SRIM, (b) LRIM. 

Finally, we need to know the value of the entropy, S q (A) at the initial integra- 
tion point A. This depends on the particular system considered, but usually 
the extreme temperature cases are known. For the Ising models (both long- 
range and short-range), Eq. (18), the limits T — > and T — > oo are: 

i. S q (T = 0) = 2™ 

ii. S q (T = oo) = 

We have implemented this Monte Carlo method using a system size iV = 34. 
Fig. 8 shows that the function that needs to be integrated in order to perform 
the T' T transformation is a smooth one. Fig. 9 compares the internal 
energy obtained by this Monte Carlo method with the exact results obtained 
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by the exact enumeration procedure showing the validity of this Monte Carlo 
scheme. 

The main disadvantage of this Monte Carlo method is that one might need 
to simulate a large range of values of T' to be able to perform accurately the 
integration needed for the T' — > T transformation. However, the use of extrap- 
olation techniques, such as the multiple histogram reweighting [24,25], which 
permit to extend the results of a simulation at a value of the temperature to a 
continuum range of temperatures, allows to reduce dramatically the number 
of simulation points needed [21]. 



5 The scaling functions 

For extensive systems, the internal energy per particle is just a function of the 
temperature, U(N,T)/N = u(T). Clearly, by definition, this scaling relation 
does not hold for non-extensive systems and there has been some recent in- 
terest in finding the correct scaling laws that apply to non-extensive systems. 
A first result was to obtain the scaling laws that follow from the application 
of Boltzmann-Gibbs statistics to a genuinely non-extensive system such as 
the LRIM in the regime a < d. The results for the internal energy, U, the 
magnetization M (defined as M = \Y^iL\Si\ and computed using Eq. (23)), 
the Helmholtz free energy F and the entropy S can be summarized by the 
following relations [26,27,18,28-30]: 



where m,u,f,s are the scaling functions. The argument justifying these scaling 
laws can be summarized as follows: the internal energy and the entropy appear 
in the definition of the free energy as F = U — TS, therefore one expects that 
U and TS should have the same behavior for large N. Since U scales as NN 
and S scales as iV one obtains that T must scale as N thus leading to the 
previous scaling ansatzs. Note that the SRIM case is recovered from the LRIM 
case in the limit a — > oo, when N — > 1 and the scaling relations, Eq. (31-34), 
become the standard ones for extensive systems. 

We present now the extensions of these scaling laws in the case that the models 
are considered under the rules of Tsallis entropy [22]. In the case q ^ 1, the 
entropy is no longer an extensive quantity (this is true both for the SRIM and 



U(N,T)=NNu(T/N) 
M(N, T) — Nm(T/N) 
F(N, T) = NN f(T/N) 
S(N,T) = Ns(T/N) 



(31) 
(32) 
(33) 
(34) 
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the LRIM). In order to generalize the argument of the previous paragraph 
giving the correct scale factor for the temperature, we derive from Eq. (3) the 
following general expression for the entropy A q (N) of a set of N independent 
particles: 

M N) = S,(N, T - oo) - I 1 + il -y - 1 (35) 

here S q (l) is the one particle entropy. In the Ising model case, one particle 
can be in any of the two states with probability 1/2, yielding: S q (l) — [1 — 
2(l/2) 9 ]/(g - 1) = [2 l - q - 1]/(1 - q). After replacing in Eq. (35), we obtain: 

A q {N) = — (36) 

Assuming that Tsallis entropy will be scaled generically with A q (N), we now 
assume that U and TS* scale in the same way as NN. Since TS/NN = 
[T A g (N) / N N][S / A q (N)] we conjecture that the temperature has to be scaled 
with N' = NN/A q (N). However, it turns out that the numerical results do not 
support this expression for the rescaling factors in the case q > 1. Therefore, 
we write the scaling relations in the following more general form: 



U q (N,T) = NNu q (T/%) (37) 

M q (N,T) = Nm q (T/%) (38) 

F q (N, T) — NNf q iT/N') (39) 

S q (N,T) = A q (N)s q (T/N' s ) (40) 

The previous argument, valid in the case q < 1, implies simply N'jj = N' s = N'. 



For consistency in the notation, we define A^(N) and A^(N) by means of 
% = NN/A q (N) and N' s = NN/A q (N), respectively. Notice that for q = 1 
it is Ai(N) oc and the scaling laws (31-34) are recovered. In order to obtain a 
good scaling description in the case g > 1 it is seen numerically that one needs 
to assume the limits A q (N) ~ 2 N ^/(q - 1) and A^(N) ~ 2 N ^/(q - 1). 
A unifying expression that reduces to the necessary ones for large iV and for 
all values of q is: 

Although we lack a satisfactory explanation for these relations, we note that 
similar scaling factors have been used previously to plot in the same scale 
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curves for the specific heat in infinite-range Ising models and non-interacting 
ideal paramagnet [31,32]. 
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Fig. 10. Internal energy U q {T) plotted in order to check the scaling relation Eq.(37) 
for q = 0.9, 1.0, 1.1 and different system sizes, (a) SRIM, (b) LRIM. 

In order to check the validity of these scaling relations, we have used the HOW 
method to simulate the one-dimensional SRIM and LRIM with a = 0.8, for 
system sizes N = 34, 100,200,400,800, 1000, and several values for the non- 
extensive parameter q e [0.1, 1.9]. We test the proposed scaling relations, Eq. 
(37-40), by plotting the scaled results in Figs. 10,11,12. One can observe in 
the Fig. 10 that, for the same value of q, the collapse of curves of different 
sizes N. This is similar to what has been observed in the q — 1 case [27]. 

It is more remarkable the fact that, with the previous choice for the scaling 
factors, all the q < 1 data collapse in a single curve. The same thing occurs for 
the q > 1 curves. Therefore, data can be described by just three universal scal- 
ing functions, corresponding to q < 1, q — 1, and q > 1 regimes respectively 
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Fig. 11. Internal energy (top graph), magnetization (middle graph) and entropy 
(lower graph) plotted in order to check the proposed scaling relations Eqs. (37,38,40) 
for the short-range Ising model (SRIM). We have use N = 1000 and the curves 
with q < 1 include q = 0.2, 0.4, 0.6, 0.8 while the curves with q > 1 include 
q = 1.2, 1.4, 1.6, 1.8. For clarity, in the entropy plot, the insert shows all the values 
of q, whereas the main plot shows only q > 1. 



(See Figs. 11,12). 
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Fig. 12. Same plot that Fig. 11 but for the long-range Ising model (LRIM). We 
plot all the same q values than in Fig. 11, although the different curves are almost 
indistinguishable with the resolution of this figure. 



One can see in Figs. 11,12 that the collapse in the entropy curves for q > 1 
is very poor. This is easily understood by noticing that the low temperature 
limit of the entropy for infinite system size is S q (T = 0) = (1 — 2 l ~ q )/(q — 1) 



26 



whereas the high temperature limit is S q (T — > oo) = 1/(5 — 1) and those 
two finite values can not be rescaled simultaneously. This is different of what 
happens for the internal energy and the magnetization for which the limits 
T — > and T — > 00 coincide for different values of q. Finally, the scaling for 
the free energy follows directly from its definition F q = U q — TS q . For q < 1 it 
is f q (x) = u q (x) — xs q (x), whereas for q > 1 and in the limit of large N, the 
scaling function is given simply by f q (x) = u q (0) — xs q (oo) = —x. 

In summary, the scaling laws given by Eqs. (37-40) work for all values of q when 
using the scaling factors given by Eqs. (41). Moreover, the scaling functions u q , 
m q and f q adopt only three different forms for each magnitude corresponding 
to q > 1, q — 1 and q < 1, both in the SRIM and LRIM cases. 



6 Thermostatistics using standard mean values 

It has been shown recently [6-8] that the use of the standard rule for the 
calculation of the mean values in Tsallis statistics provides also a valid ther- 
modynamical formalism. By "standard" rule we mean the use of the first 
option for the averages in which u(pi) = Pi is used in (4). Moreover, it has 
been argued [4,33] that the results of using this first option coincide with the 
results of the third option (the one used up to here in this paper) with a trivial 
change in the parameter q — > 1/q. In this section we show that it is possible 
indeed to map the results of one option into the results of the other, although 
the relation between them implies, besides the previous change in the param- 
eter q, a non-trivial mapping for the temperature. Numerical results using the 
techniques developed in the previous sections, will allow us to plot the relation 
between the temperatures of the two options. 

For the sake of clarity in the exposition we will use the subindexes "1" and 
"3" to denote the results one obtains in each option. Hence, the first option 
seeks the maximization of 

S 1 (q) = 1 " E f (42) 
q - 1 

subject to the canonical ensemble constrains: Pi > 0, J2iPi — 1, Hi^iPi = U. 
The third option, on the other hand, seeks the maximization of 

S 3 (q) = ^ > (43) 

q - 1 

subject to the same constrains. Of course, in the third option, the probabilities 
Pi should be interpreted as "escort" probabilities, but this interpretation has 
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no practical consequence whatsoever in the calculation of the averages. The 
key point now is that both entropies are related by: 

S 1 (l/q) = G q [S 3 (q)} (44) 




Fig. 13. Transformation function between the two entropic forms S 3 (q) and Si(l/q) 
as indicated in Eq. (44). 

Where G q (x) = — (1 + (1— q)x)~ l l q ] is a monotonically increasing function 
of x. This function satisfies the property: G~ 1 {x) = Gi/ q (x), see Fig. (13). 
Hence, the same set of probabilities {p^ that maximize S 3 (q) for a given 
value of U will maximize Si(l/q), for the same value for U. However, the 
fact that the probabilities coincide in both options does not mean that the 
averages computed using these probabilities coincide when they are plotted 
as a function of the temperature because it turns out that there is a non- 
trivial relation between the temperatures of both options. Let us denote by 
Ti and T 3 , respectively, the temperatures of the 1st and 3rd options. They 
can be defined as the (inverse of the) Lagrange multiplier needed to satisfy 
the constraint of fixed mean energy or, alternatively, they have been shown to 
satisfy the relations [7]: 

i/r lh ) = ^ 1/T,(,) = (45) 

Using (44) we find the desired relation between the two temperatures: 

Ts(q)=T 1 (l/q)/G' 1/q [S 1 (l/q)} (46) 

G' (x) is the derivative of G q (x). After substitution of Eq.(42), we find: 

TM=T l {l/q)(£pl' q ) q+1 (47) 
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Therefore, it is true that the results of the third option at the value q of 
the parameter can be obtained from those of the first one at the value 1/q. 
However, the mapping requires a non-trivial rescaling of the temperature, as 
given by Eq.(47). Let us recall again that only the dependence with T does not 
vary when changing the zero of energies and, hence, can be the only physically 
relevant one. 
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Fig. 14. The same plot as in Fig. (7) but using lineal mean values instead non-lineal 
ones for the values of N indicated and q = 0.8, 0.9, 1.0, 1.1, 1.2. (a) SRIM, (b) LRIM. 



In order to give an alternative explanation of the relation between the tem- 
peratures of both options, let us write down the solutions for the probabilities 
using the f3' parameter. For the third option, the solution is read directly from 
Eq.(17) that we rewrite using the notation of this section: 

0, l-(l-g)fe<0 

Vl = { n-a-w^, , otherwise (48) 
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Fig. 15. Transformation between temperatures T\ and T3 using Eq. (47) for N = 100 
and the indicated q values, (a) SRIM, (b) LRIM. 

where the parameter /3' 3 = I/T3 is related to the temperature T 3 by Eq.(16) 
which reads: 



T'M = T 3 (q)[l + (1 - q)S 3 (q)} + (1 - g)f/ 



(49) 



For the first option, it is possible to write the solution in the following form: 

0, 



Pi = < 



[i-(i-i/ g )fl e< ]g- 



1 - (1 - < 

otherwise 



(50) 



where the parameter f3[ = l/T[ is related to the temperature 7\ by [34] 

T[(g) = T!(g)[l + (1 - q)S 1 (q)] + (1 - 1/ ? )C/ (51) 
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Now, it is clear by comparing Eqs. (48) and (50) that the probabilities pi of the 
third option computed at q are equal to the probabilities Pi of the first option 
computed at 1/q provided that we choose the same values for the primed 
temperatures, T^(q) = T[(l/q). After substitution of (49) and (51) and using 
(43), (42) we recover Eq.(47) 

It is straightforward now to use the number of states fl(Ek) obtained using 
the HOW method to compute Eqs. (50,51,42) by replacing the sums over the 
configurations to sums over energy levels weighted by VL(E k ). In this way, 
we can perform the necessary averages implied in the first option as well as 
the temperature transformation factor needed in Eq. (47). In Fig. 14 we plot 
the internal energy U q as a function of the temperature using the standard 
averages of the first option. The most noticeable difference with the results 
of the third option, see figure (5) is that it is not necessary now to use the 
Maxwell construction because there are no loops with the temperature. In 
Fig. 15 we plot T 3 (q) vs. T^l/q) in the LRIM and SRIM cases. Using these 
two results, it is possible to obtain the averages within the third option as 
a function of T3. Of course, the results agree perfectly with those shown in 
Fig. 7. It is possible also to obtain from Eqs. (37-40) the scaling relations valid 
when the standard calculation of mean values is used for the calculation of 
thermodynamic quantitities. 



7 Microcanonical Ensemble 



As mentioned in the introduction, the third option can be formulated by using 
the entropic form Eq. (7) plus the standard rule for the calculation of mean 
values, Eq. (5) or, alternatively, by using the original entropic form Eq. (2), 
but with a mean value definition: 

(0) q = T,Oi^hi- (52) 



These two points of view are completely equivalent. The first option, as ex- 
plained in the previous section, uses also the original entropic form but with 
standard mean values. We will consider in this section Tsallis original entropic 
form S(q) in the context of the microcanonical formalism. The aim is to be 
able to derive the internal energy without any a priori assumption about the 
definition of averages. In the microcanonical ensemble we consider the maxi- 
mization problem for the original entropic form S q given by Eq.(2) with the 
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constrain of given energy E. The solution is the equiprobability, 

n(E)~\ Si = E 



Pi = s 



0, 



otherwise 



(53) 



Where Q(E) is the number of configurations with energy E. The entropy as 
a function of the energy is: 



Q{E) l - q - 1 
1^ 



(54) 




200 - 



Fig. 16. Plot of the internal energy as a function of the temperature for N = 100. 
By points we show the results of using the microcanonical ensemble described in 
the text, and by lines the results obtained from the canonical ensemble (see Fig. 7) 
by using the non standard mean values of the third option, (a) SRIM, (b) LRIM. 
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The temperature is defined by the thermodynamic relation Eq. (13), ^ = ^ 
or 



T K J dE 



(55) 



Inverting this relation, we obtain the energy as a function of the temperature, 
E(T). In general, this relation needs to be inverted numerically. In terms of 
the scaling function <j)(x) defined in (25) we have: 



T = 



e (q-l)N4>(E/NN) 

f(E/NN) 



(56) 
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Fig. 17. Plot of the energy fluctuations in the canonical ensemble 
cr(7i)/NN = yj (TC 2 ) — (7i)' 2 /NN as a function of the scaled temperature for the 
indicated values of the sizes N and the q parameters, (a) SRIM, (b) LRIM. 

where (f>(x) is known exactly for the SRIM and can be evaluated numerically 
using the HOW method for the LRIM (from the plot in the insert of Fig. 5b). 
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Results are shown in Fig 16 where we plot the internal energy coming from the 
application of the microcanonical ensemble to both the SRIM and the LRIM. 
In the same figure we have also plotted the energy coming from the canonical 
ensemble using non standard mean values (third option). We can see that both 
approaches coincide for q < 1, and differ for q > 1 in some temperature range. 
The ultimate reason for not having equivalence between the two ensembles 
is that fluctuations of the energy in the canonical ensemble can not be ne- 
glected. We have checked that this is indeed the case by computing the energy 



fluctuations o"(7i) = \j (H 2 ) — (H) 2 as a function of the system size. In the 

Fig. (17) we see that fluctuations, normalized by the scale of energy, NN, do 
not decay to zero for increasing N in the range of temperatures for which the 
microcanonical and canonical ensemble do not agree. For q < 1 fluctuations 
do decay to zero with the system size in all the temperature range. 

If we compare the microcanonical and the canonical ensemble using the stan- 
dard mean values (the first option studied in the last section) we observe the 
coincidence of both ensembles for q > 1, and disagreement (in a given temper- 
ature range) for q < 1. This turns out to be also consistent with non-vanishing 
energy fluctuations in the appropriate range. This is the expected result be- 
cause of the mapping q — ■> 1/q implied in going from the third option to the 
first one. 



8 Conclusions 

In this paper, we have given details of two methods that can be used to 
perform numerical simulations for many particle systems that are governed 
by generalized statistics, such as the Tsallis one. The first method extends 
the Histogram by Overlapping Windows method, devised originally for short 
range Hamiltonians, to systems with very-long range interactions. The sec- 
ond method, devised specifically for the Tsallis thermostatistics, uses a typical 
Metropolis Monte Carlo updating scheme combined with a numerical integra- 
tion. We have emphasized the need of using the right temperature definition 
if averages are to be independent of the zero of energy. We have applied our 
methods to the case of the Ising model with either short range (SRIM) or 
long range (LRIM) interactions. The latter case corresponds to a situation 
genuinely non-extensive in which the energy levels scale as N x with x > 1, 
N being the number of variables. We have compared the methods with some 
exact results available in the case of one-dimensional short range Ising model 
of arbitrary size and the long-range model for small system sizes. 

We have shown that the internal energy, entropy, Helmholtz free-energy and 
magnetization follow non trivial scaling laws with the temperature T and the 
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number of variables N. We have justified these scaling laws by some heuristic 
arguments that, however, fail to reproduce the observed behavior for q > 1. 
These scaling laws for q ^ 1 are non-extensive in the sense that the different 
thermodynamic potentials have to be scaled with a factor that depends in a 
non-trivial, i.e. non linear, way of N. The scaling laws hold for both the LRIM 
and the SRIM (with different scaling functions in each case), independently 
of the fact that the systems are genuinely non-extensive or extensive. This 
shows that the non-extensivity arises mainly because of the application of the 
Tsallis statistics. 

We have discussed the differences between the use of standard (first option) 
and non-standard (third option) mean value definitions for the Tsallis Ther- 
mostatistics formalism. We show that, although the results of both definitions 
can be mapped onto each other by using the q — > 1/q transformations, this 
mapping requires as well a non-trivial change in the temperature. Finally, we 
have shown that the use of the microcanonical ensemble coincides with the 
results of the canonical ensemble in the third option only for q < 1. We in- 
terpret this result as the non-vanishing energy fluctuations that occur in the 
corresponding case. 

An obvious extension of the results presented here is to consider the Ising 
models in spatial dimension greater than one, [35]. The HOW method can be 
extended in any dimension for short range and long range interactions. For 
the SRIM in d = 2 the exact results [36] should be used. 

We remark that the present work concerns equilibrium systems and there is no 
time dependence in our simulations. However, it has been recently conjectured 
[3] that Tsallis statistics appears in some non-equilibrium systems such as the 
relaxation of the non-neutral plasma experiments in [37]. To study these non- 
equilibrium systems within the Tsallis statistics formalism, it would be more 
appropriate to use Molecular Dynamics (MD) methods in which the evolution 
equations are solved as a function of time. We are currently working on a MD 
simulation valid for a Lennard- Jones system within the Tsallis statistics. This 
MD method uses a Kusnezov, Bulgac and Bauer thermostat [38,39] where ad- 
ditionally the actual temperature has to be calculated using a relation similar 
to the Eq. (30). 
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